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We provide a comprehensive study of the cosmic-ray muon flux and induced activity as a func¬ 
tion of overburden along with a convenient parameterization of the salient fluxes and differential 
distributions for a suite of underground laboratories ranging in depth from ~1 to 8 km.w.e.. Partic¬ 
ular attention is given to the muon-induced fast neutron activity for the underground sites and we 
develop a Depth-Sensitivity-Relation to characterize the effect of such background in experiments 
searching for WIMP dark matter and neutrinoless double beta decay. 
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I. INTRODUCTION 

Underground laboratories provide the overburden nec¬ 
essary for experiments sensitive to cosmic-ray muons and 
their progenies. Muons traversing a detector and its sur¬ 
rounding material that miss an external veto serve as a 
background themselves and secondary backgrounds are 
induced in the production of fast neutrons and cosmo¬ 
gonic radioactivity. In this study we have focused on the 
muon-induced fast neutron background as a function of 
depth and the implications for rare event searches for 
neutrinoless double beta decay and WIMP dark matter. 
One of our main goals is to develop a Depth-Sensitivity- 
Relation (DSR) in terms of the total muon and muon- 
induced neutron flux and to put this into the context of 
existing underground laboratories covering a wide range 
of overburden. 

In Section II we review the experimental data avail¬ 
able for differential muon fluxes and provide a definition 
of depth in terms of the total muon flux that removes 
some confusion regarding the equivalent depth of an un¬ 
derground site situated under a mountain versus one with 
flat overburden. The muon fluxes and differential distri¬ 
butions are parameterized and used as imut in Section 
III to generate, via FLUKA simulations [y, the produc¬ 
tion rate for fast neutrons. The total neutron flux and 
salient distributions are compared with the available ex¬ 
perimental data and we provide some convenient parame- 
terizations that can be used as input for detector-specific 
simulations at a given underground site. We quantify the 
agreement between FLUKA simulation and experimen¬ 
tal data and provide an explanation for the discrepancy 
between neutron flux and energy spectra as measured in 
the LVD detector. Muon-induced cosmogenic radioactiv¬ 
ity is discussed in terms of depth and the average muon 
energy in Section IV. In Section V we apply our results 
to a generic study of germanium-based experiments in 
search of neutrinoless double beta decay and WIMP dark 
matter and demonstrate the utility of the DSR in pro¬ 
jecting the sensitivity and depth requirements of such 
experiments. We conclude with a summary of the results 


and an outline of new studies underway. 


II. DEPTH-INTENSITY-RELATION AND 
DISTRIBUTIONS FOR COSMIC-RAY MUONS 

A. Through-Going Muon Intensity 

1. Differential Muon Intensity versus Slant-Depth 

The cosmic-ray muon flux in the atmosphere, under¬ 
ground, and underwater has been a subject of study for 
more than five decades [^. Experimental data on the 
differential muon intensity (in units of cm~'^s~^sr~^) are 
shown in Fig. ^ as a function of slant-depth measured in 
units of kilometers of water equivalent (km.w.e.), where 
1000 hg/cm^ = 10® g/cm^ = 1 km.w.e.. 

Groom et. al. proposed a model [it'll to fit the experi¬ 
mental data to a Depth-Intensity-Relation (DIR), appro¬ 
priate for the range (1 - 10 km.w.e.): 

I{h) = -b /2e(-^/^=)), (1) 

where I(/i) is the differential muon intensity correspond¬ 
ing to the slant-depth, h. 

Using the experimental data in Fig. ^ we deter¬ 
mine the free parameters of equation o as: Ii = 
(8.60±0.53)x 10-® sec-Um-2sr-i, I 2 = (0.44±0.06)x 
10“® sec“^cm“^sr“^, Ai = 0.45±0.01 km.w.e., A 2 = 
0.87±0.02 km.w.e.. The relative deviation between the 
data and our fit is shown in Fig. (21 indicating that the 
parameterization reproduces the experimental data rea¬ 
sonably well and with an overall accuracy of about 5%. 


2. The Total Muon Flux with Flat Overburden 

For an underground laboratory with flat overburden it 
is straightforward to calculate the total muon intensity 
arriving below the surface at a vertical depth, /iq. In the 
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FIG. 1: Measurements of the differential muon flux as 
a function of slant depth from Castagnoli Q. Barrett H, 
Miyake 0 , WIPP 0 , Soudan 0 , Kamioka Boulby 0 , 
Gran Sasso mill, Frejus HI and Sudbury fl^ . Note that 
the measurements for Kamioka 0 and Sudbury are re¬ 
ported as the number of muons per day. We calculate the 
effective detector acceptance for these two measurements in 
order to obtain the muon flux. The solid curve is our global 
fit function described by equation Q- 



FIG. 2: The relative deviation between the global fit func¬ 
tion and the measured data on the differential muon flux from 
Gastagnoli 0 , Barrett 0 , Miyake 0 . WIPP 0 . S oudan 0 , 
Kamioka 0 , Boulby 0 , Gran Sasso K(l||. Freiusll3 and Sud¬ 
bury m. The horizontal lines indicate the root-mean-square 
deviation amongst the residuals. 

flat earth approximation, the through-going muon inten¬ 
sity (lift,) for a specific slant-depth, h, in the direction of 
zenith angle, 9, reads: 

Ith{h,9)=I{h)G{h,9), (2) 

where G{h,9) = sec(0), h = /iosec(0), and 1(h) is the DIR 
expressed in equation O- Equation o now becomes 

Ith(h,e) = 

( 3 ) 


Integration over the upper hemisphere using equation 
© then provides the total muon intensity for an under¬ 
ground site with flat overburden positioned at a vertical 
depth /lo¬ 
using the experimental data for the total muon flux 
and knowledge of the vertical depth for a set of under¬ 
ground sites with flat overburden (WIPP 0 , Soudan 00 , 
Boulby 0 and Sudbury 0 ) we can now define a fit- 
function which is similar to the differential muon inten¬ 
sity function (equation (^l): 

I^(ho) = 67.97 X -b 2.071 x 10-®e^, (4) 

where /iq is the vertical depth in km.w.e. and is in 

units of cm“^s“^, appropriate in the flat-earth approxi¬ 
mation. 

3. The Total Muon Flux in Case of Mountain Overburden 

In the case that a laboratory is situated underneath 
a mountain, additional information regarding the moun¬ 
tain shape or elevation map, h(0, (/>), is required to deter¬ 
mine the total muon flux: 

Itot = J Sin(9)d9 J dct)I(h(e,cj)))G(h, 9), (5) 

where G(h,9) = sec(9) and Itot is the total muon flux 
obtained after integrating over the mountain shape and 
using the DIR defined in equation ©• 

As an example, we have computed the total muon flux 
at the Gran Sasso Laboratory using the detailed infor¬ 
mation provided by the MACRO collaboration 0 on 
the mountain shape and their measurements of the dif¬ 
ferential muon flux (see Fig. 0. We And a total muon 
intensity of (2.58±0.3) x 10“®cr7i“^sec“^, which is consis¬ 
tent within about 20% to that obtained in Refs. 00 . 
If this intensity is now entered into the left hand side of 
equation 0, we can now solve for the equivalent vertical 
depth relative to a flat overburden for the Gran Sasso 
Laboratory and And it to be 3.1 ± 0.2 km.w.e.. 

This depth should not be confused with the average 
depth that would be deduced simply by integrating over 
the depth profile of the mountain: 

<h>= J sin(9)d9 J d(j)h{9, cj)), ( 6 ) 

which yields 3.65 km.w.e.. A similar approach can be 
taken with information available from the Frejus Col¬ 
laboration 0 . We And a total muon intensity of 
(4.83±0.5) X 10“®cm“^sec“^ corresponding to an equiva¬ 
lent flat overburden of 4.2 ± 0.2 km.w.e. and an average 
depth of 5 km.w.e.. Our calculation is consistent with 
the Frejus Collaboration’s result within 12%. We note 
that the equivalent “flat-overburden” depth defined by 
the experimental measure of the total muon flux is ~15- 
20% lower than that often quoted for Gran Sasso and 
Frejus based on the average physical depth. 
















J^. Definition of Depth and Total Muon Flux for 
Underground Sites 
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The data on the total muon intensity at the various 
underground sites is summarized in Table d and Fig. 13 
We use equation 0 to calculate the total muon flux 
for Homestake (flat-overburden) at the depth 4.3 ± 0.2 
km.w.e. ^ 3 . The relative difference between the data 
and our model (equation Q)) is shown in Fig. ^ where 
the uncertainties reflect the experimental uncertainties 
in Table ^ In order to circumvent the misuse of vertical 
muon intensity in comparing sites with flat overburden 
to those under mountains, we define the equivalent depth 
relative to a flat overburden by the experimental mea¬ 
surements of the total muon intensity. This definition 
and these intensities are used hereafter. 



FIG. 3: The total muon flux measured for the various un¬ 
derground sites summarized in Table Q] as a function of the 
equivalent vertical depth relative to a flat overburden. The 
smooth curve is our global fit function to those data taken 
from sites with flat overburden (equation 0). 


TABLE I: Summary of the total muon flux measured at the 
underground sites and the equivalent vertical depth relative 
to a flat overburden. 


Site 

Total flux 
cm“^sec~^ 

Depth 

km.w.e. 

WIPP 

(4.77T0.09) X 10~' 

1.585T0.011 

Soudan 

(2.o±o.2) X 10 "^ yj 

1.95T0.15 

Kamioka 

(1.58±0.21) X 10“'^ [2] 

2.05±0.15t 

Boulby 

(4.09±0.15) X 10~® 

2.805T0.015 

Gran Sasso 

(2.58±0.3) X 10“®[this work] 

3.1±0.2t 


(2.78±0.2) X 10"® [16] 

3.05±0.2t 


(3.22±0.2) X 10"® ] 17 ^ 

2.96±0.2t 

Frejus 

(5.47±0.1) X 10"® [14] 

4.15±0.2t 


(4.83 ±0.5) X 10"® [this work] 

4.2±0.2t 

Homestake 

(4.4±0.1 X 10"®) [this work] 

4.3T0.2 

Sudbury 

(3.77T0.41) X 10"^° [12] 

6.011T0.1 


^ Equivalent vertical depth with a flat overburden 
determined by the measured total muon flux. 
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FIG. 4: The relative deviation between data on the total 
muon flux and our global fit function. The horizontal lines 
indicate the root-mean-square deviation amongst the resid¬ 
uals based upon the experimental uncertainties in the mea¬ 
surements. 


B. Stopping Muon Intensity 

Stopping-muons are also a source of background. For 
example, capture on a nucleus produces neutrons 
and radioactive isotopes. The total stopping-muon rate 
has contributions from cosmic-ray muons coming to 
the end of their range, secondary muons generated lo¬ 
cally through interactions of the primary muons (due to 
virtual-photo interactions with nuclei), and local muon 
production by real photons (TTg-decay in electromagnetic 
showers). It is customary to quote results in terms of the 
ratio, R, of stopping muons to through-going muons. A 
detailed calculation is provided by Cassiday et al. [l^ . 
The total ratio, R(/i), of stopping-muons to through- 
going muons (vertical direction) at different depths can 
be parameterized as [T^ 

AEe'^/^ 

where 7 ^ = 3.77 for > 1000 GeV 1^. ^ = 2.5 km.w.e., 
AE ^ ah, a = 0.268 GeV/km.w.e. |2ll for E^ > 1000 
GeV I 23 , h is the depth of an underground laboratory, 
and €fi = 618 GeV [^. For large depths, as can be seen 
in Fig. 13 this ratio is less than 0.5% and is hereafter 
neglected for the underground sites considered in this 
study. 


C. Muon Energy Spectrum and Angular 
Distribution 

In addition to the total muon intensity arriving at a 
given underground site, we require knowledge of the dif¬ 
ferential energy and angular distributions in order to gen¬ 
erate the muon-induced activity within a particular ex¬ 
perimental cavern. The energy spectrum is discussed in 
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FIG. 5: The ratio of stopping-muons to through-going muons, 
relative to the vertical direction, as a function of depth. 


Refs. [Hil: 


dN 

dEf_, 




( 8 ) 


where A is a normalization constant with respect to the 
differential muon intensity at a given depth and is 
the muon energy after crossing the rock slant depth h 
(km.w.e.). Fig.^lshows the local muon energy spectrum 
for the various underground laboratories under consider¬ 
ation using the parameters b = 0.4/km.w.e., 7^ = 3.77 
and = 693 GeV |^. Fig. [ 7 | shows the local angular 
distribution for the same sites where we assume a secjO) 
distribution, valid for depths in excess of 1.5 km.w.e. |24| . 
Note that the overall angular distribution of muons at the 
surface is proportional to cos‘^{9) with an average muon 
energy of about 4 GeV . 
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FIG. 6: The muon energy spectrum local to the various un¬ 
derground sites calculated using equation ®. The areas un¬ 
der the curves are normalized to the vertical muon intensity 
for comparison purposes. 

From equation 0 , the average muon energy at depth 



FIG. 7: The muon angular distribution local to the various 
underground sites based on equation 0. All curves have 
been normalized to the total muon intensity for comparison 
purposes. 


h is given by: 


< >= 


7m - 2 


(9) 


The parameters e^, b and 7,,, in equation have been 
studied by several authors for standard rock 

(A = 22, Z = 11, p = 2.65 g cm“^). Uncertainty in these 
parameters are due to uncertainties in the muon energy 
spectrum in the atmosphere, details of muon energy loss 
in the media, and the local rock density and composi¬ 
tion. Table HD summarizes the average muon energy for 
the various sites where we have used two different sets of 
parameters provided by Lipari et al. ( b = 0.383/km.w.e., 
7^ = 3.7 an d e,, = 618 GeV [^) and Groom et al. (b = 
0.4/km.w.e. [H, 7^ = 3.77 and = 693 GeV HI). The 
measured average single muon energy at Gran Sasso HI 
is 270±3(stat) ± 18(syst) GeV which has an uncertainty 
of 6.8%. The predicted values using both sets of param¬ 
eters agree with the measured value within the measured 
uncertainty. 


TABLE II: Single muon average energies for the various un¬ 
derground sites. 


Site 

Lipari et al. 

Groom et al. Measured value 

WIPP 

165 GeV 

184 GeV 

Soudan 

191 GeV 

212 GeV 

Kamioka 

198 GeV 

219 GeV 

Boulby 

239 GeV 

264 GeV 

Gran Sasso 253 GeV 

278 GeV 270T18 GeV [2^ 

Sudbury 

327 GeV 

356 GeV 


III. MUON-INDUCED NEUTRONS 

We distinguish two classes of fast neutrons, namely 
neutrons produced by muons traversing the detector it- 
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self, and neutrons created in the external rock by muons 
missing the veto detector. The former can be tagged ef¬ 
fectively in an external veto with sufficient efficiency sur¬ 
rounding a central detector. The latter are more difficult 
to shield or tag in coincidence with the primary muon 
owing to the hard energy spectrum and long propagation 
range. Thus, we focus here on the fast neutrons produced 
in the external rock and quantify the production rate as 
a function of depth. 

The production of fast neutrons depends strongly on 
the depth and composition of an underground site. Gen¬ 
erally speaking, the neutron production rate at large 
depths due to muons is two to three orders of magnitude 
smaller than that of neutrons arising from local radioac¬ 
tivity through (a,n) reactions. Nonetheless, the latter 
process is common to any underground experiment and 
the low energy neutrons (typically <8 MeV) produced via 
(a,n) reactions are relatively easy to shield (see Section 
V.D). The muon-induced neutrons, on the other hand, 
have a very hard energy spectrum (extending to several 
GeV) and can penetrate to significant depth both in the 
surrounding rock and detector shielding materials. In 
this section we exploit the total muon fluxes and dis¬ 
tributions developed in the previous section as input to 
FLUKA simulations to study the muon-induced neutron 
flux, energy spectrum, angular distribution, multiplicity, 
and lateral distribution in the underground laboratories. 

Table lini exhibits the rock composition for the six sites 
under consideration. The rock composition of the WIPP 
site is mainly NaCl [^, whereas for the other four sites, 
the average atomic weight and average atomic number 
are calculated based on the known local rock composi¬ 
tion nimni. For lack of additional information we 
assume standard rock for Kamioka. Note that Boulb y is 
a salt mine but the rock composition provided in Ref. 
is very similar to the standard rock. 


TABLE III: Average matter properties of the various under¬ 
ground sites. 


Site 

A 

V 

A 

84 

V 

< Z >/< A > g/cm^ 

WIPP 

30.0 

14.64 

0.488 

2.3 

Soudan 

24.47 

12.15 

0.497 

2.8 

Kamioka 

22.0 

11.0 

0.5 

2.65 

Boulby 

23.6 

11.7 

0.496 

2.7 

Gran Sasso 22.87 

11.41 

0.499 

2.71 

Sudbury 

24.77 

12.15 

0.491 

2.894 


The rock thickness employed in the simulation is 20 m 
X 20 m X 20 m. The laboratory cavern of size 6 m x 
6 m X 6 m was placed inside the rock region at a depth 
of 7 m from each side of the cube. This cube ensures 
equilibrium between neutron and muon fluxes, hence the 
ratio of neutron to muon fluxes is constant. 


A. Comparison Between Data and Simulation 

Table Hvl lists the mean neutron pro duction rates from 
seven measurements using liq¬ 

uid scintillator covering a significant range in depth and 
mean muon energy. We provide a global fit function to 
the data as a function of mean muon energy (see Fig. El) 
and compare this to the Monte Carlo calculations per¬ 
formed in Ref. (C 10 H 22 ), Ref. (C 10 H 20 ) and our 
FLUKA simulation (C 10 H 20 ). For experiments which did 
not provide the mean muon energy, we use the experi¬ 
mental depth and the muon energy loss rate [U to esti¬ 
mate the mean muon energy. 


TABLE IV: Measured neutron production rates. 


Measurements 

Depth 

km.w.e. 

0 A 

V 

< n > 

n/{iigcm~^) 

Hertenberger 1291 

0.02 

13 

(2T0.7) X lO”'" 

Bezrukov 

0.025 

14.7 

(4.7T0.5) X 10“® 

Boehm 

0.032 

16.5 

(3.6±0.31) X 10“® 

Bezrukov 

0.316 

55 

(1.21T0.12) X lO"'^ 

Enikeev 

0.75 

120 

(2.15±0.15) X 10"'‘ 

the LVD data [3^ 

3.1 

270 

(1.5±0.4) X 10“"* 

Aglietta (the LSD) [33] 

5.0 

346 

(5.3±1.1) X lO”* 


Note that the LVD (see Fig. El and Table IIVI) result 
obtained at Gran Sasso deviates significantly from the 
global-fit curve and simulations. Not only is the mea¬ 
sured flux apparently low, the differential energy spec¬ 
trum of fast neutrons measured in the same experiment is 
also inconsistent with simulation. Kudryavtsev et al. 
suggested that the quenchin g of proton-recoil energy in 
the liquid scintillator of LVD EJ] is a natural explanation 
for the discrepancy between simulation and the measured 
energy spectrum. We propose that this same effect is also 
responsible for the discrepancy in the measured neutron 
flux. 

Following directly from Ref. describing the LVD 
analysis for neutrons, “a high-energy-threshold (HET) 
trigger is set at 4-5 MeV. During the 1 ms time period 
following an HET trigger, a low-energy-threshold (LET) 
is enabled for counters belonging to the same quarter 
of the tower which allows the detection of the 2.2 MeV 
photons from neutron capture by protons. Each neutron 
ideally should generate two pulses: the first pulse above 
the HET is due to the recoil protons from n-p elastic 
scattering (its amplitude is proportional to and even close 
to the neutron energy); the second pulse, above the LET 
in the time gate of about 1 ms is due to the 2.2 MeV 
gamma from neutron capture by a proton. The sequence 
of two pulses (one above HET and one above LET) was 
the signature for neutron detection. The energy of the 
first pulse (above HET) was measured and attributed to 
the neutron energy.” 

The authors go on to “Note that really this is not the 
neutron energy but the energy transfered to protons in 
the scintillator and measured by the counter”, however. 
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they do not correct the visible energy for quenching ef¬ 
fects to yield the true proton-recoil energy. Due to the 
finite HET, this also means that the total number of neu¬ 
trons counted is also underestimated. Quenching of pro¬ 
tons in scintillator was measured by Ref. and a factor 
of 2.15 is expected for 4 MeV energy of proton recoil as 
shown in Fig.^ If we correct the LVD results by taking 
this quenching factor and energy threshold into account 
the total neutron production rate we obtain is 4.5x10“^ 
n/(/i g cm~‘^), an increase of about a factor of three from 
the published value and consistent with our global 
fit curve. 

One should clarify that such a correction does not ap¬ 
ply to the LSD data 13.111 taken from an experiment sim¬ 
ilar to LVD and operated at Mont Blanc Laboratory. In 
LSD, no attempt was made to measure the energy of the 
muon-induced neutrons, however, neutrons were counted 
by demanding a HET (25-30 MeV) produced by muons 
with a track length of at least 15 cm in the liquid scintilla¬ 
tor. Neutrons were then tagged in coincidence with the 
delayed capture gamma-ray. Consequently, apart from 
minor corrections owing to those initial muons produc¬ 
ing coincident neutrons that miss the muon trigger, there 
should be no significant threshold correction associated 
with neutron counting in LSD. 



FIG. 8: The neutron production rate in liquid scintillator 
versus the mean muon energy. Data points with uncertain¬ 
ties are experimental measurements from Hertenberger [23, 
Boehm l3fl . Bezrukov 1^ . Enikeev l3^ . the LVD datal^ 
and Aglietta . The solid curve is our global fit to the data 
after correcting the LVD data point for quenching effects de¬ 
scribed in the text. Our global fit curve describes the data 
well but the FLUKA simulations tend to underestimate the 
neutron production rate by about 35%. 

As can be seen in Fig’s |S1 and d the data are well 
described by a simple power law model suggested by 
Refs. [ 3 ^ 1^ [33 and our FLUKA simulation. The 
FLUKA simulations, however, underestimate the data 
(and the simple power law fit to this data) by about 30 
percent. It is natural to attribute this to the virtual 
photonuclear cross section which is not well known for 
high energy cosmic-ray muon interactions with nuclei. 



FIG. 9: The quenching function relating the effective or visi¬ 
ble energy in liquid scintillator versus the kinetic energy im¬ 
parted to a recoiling proton induced by neutron scattering 
(adopted from Ref. I^ i. 



10^ 

Mean Muon Energy (GeV) 


FIG. 10: The relative deviation between our global fit func¬ 
tion and the measured neutron production rate as a function 
of the mean muon energy. The uncertainties represent the 
experimental uncertainties on the measured data points and 
the horizontal lines indicate the root-mean-square deviation 
on the residuals. 


Nonetheless, the integrated cross section of virtual pho¬ 
tonuclear interactions of muons measured by MACRO 
and ATLAS show agreement with the Bezrukov- 
Bugaev model used in FLUKA, though the accuracy 
of the prediction is limited by the lack of data for muon- 
induced interactions in materials of medium density and 
composition. 

We suggest that it is possible that the neutron multi¬ 
plicity in the muon-induced nuclear cascades and electro¬ 
magnetic cascades is responsible for this difference. The 
experimental results from Refs. show a higher 

neutron multiplicity than that predicted by FLUKA 
and we propose a neutron multiplicity correction func¬ 
tion to correct the neutron production rate. This func¬ 
tion is obtained by extrapolating the variation in neu- 
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tron multiplicity as a function of muon energy between 
the proposed parameterization based on the measure¬ 
ments 1^ and the FLUKA simulation. 


Md - Mmc 

Md 




( 10 ) 


where Md is the measured neutron multiplicity, Mmc is 
the simulated neutron multiplicity in FLUKA and is 
the muon energy in GeV. After correcting the neutron 
multiplicity in the FLUKA simulation, good agreement 
is found between the data and the simulation as can be 
seen in Fig. m 



FIG. 11: The muon-induced neutron production rate versus 
the mean muon energy after correcting the neutron multiplic¬ 
ity in the FLUKA simulation. 

Further improvement might be gained with minor 
modifications in the inclusion of deep inelastic scatter¬ 
ing of muons on nucleons. More generally, it is desir¬ 
able to have more data for high energy muon interac¬ 
tions in the appropriate materials in order to more ac¬ 
curately tune the simulations relevant to neutron pro¬ 
duction deep underground. Nonetheless, whether we use 
our global fit function to the measured data or rely upon 
the multiplicity-corrected FLUKA simulation, the muon- 
induced neutron yield reproduces the data within an ac¬ 
curacy of about 15%. 


B. Media Dependence of Neutron Production Rate 

The muon-induced production rate for neutrons de¬ 
pends critically on knowledge of the chemical compo¬ 
sition and density of the medium through which the 
muons interact. We have studied this dependence us¬ 
ing the FLUKA simulation specific to Gran Sasso in or¬ 
der to compare directly with Ref. The dependence 
on atomic weight is shown in Fig. 112 where the general 
trend is well described by a power law, consistent with 
Ref. using slightly different fitting parameters. 

<n>= 4.54 X 10“®A° ®in/(Ai g cm-"^). (11) 



FIG. 12: Simulation of the muon-induced neutron production 
rate versus the atomic weight of the medium. 


The contribution to the neutron production rate from 
electromagnetic showers becomes more important for a 
heavy target, since the cross-section of an electromag¬ 
netic muon interaction is proportional to T? jA. Fig. 
shows this dependence where, again, the general trend 
can be described using a power law: 

<n>= 1.27 X 10-4(:^)O-92j^/(^ g cto-2). (12) 



FIG. 13: Simulation of the muon-induced neutron production 
rate versus T? jA of the medium. 


C. Neutron Fluxes and Differential Spectra at 
Underground Sites 

1. Neutron Flux at Rock/Cavem Boundary 

The muon-induced neutron flux emerging from the 
rock into the cavern has been estimated for the various 
underground sites considered in this work. We derive 
the neutron flux utilizing the FLUKA simulation with 









the corrected neutron multiplicity (equation llOll 'l and 
the muon fluxes and distributions outlined in Section 11. 
The neutron flux {4>n) as a function of depth is shown 
in Fig. m where we have included a fit function of the 
following form: 

= (13) 

where ho is the equivalent vertical depth (in km.w.e.) 
relative to a flat overburden. The fit parameters are Pq = 
(4.0 ± l.l)xl0“^ cm“^s“^ and Pi = 0.86 ± 0.05 km.w.e.. 



FIG. 14: The total muon-induced neutron flux deduced for 
the various underground sites displayed. Uncertainties on 
each point reflect those added in quadrature from uncertain¬ 
ties in knowledge of the absolute muon fluxes and neutron 
production rates based upon our simulations constrained by 
the available experimental data. 

In Table El we summarize the neutron flux at the 
rock/cavern boundary for the various sites considered 
and note that we have not included the effect of neutrons 
that emerge from one surface and back-scatter back into 
the cavity. The results are in good agreement with the 
existing simulation results for Gran Sasso ^3- If 
simulation results for Boulby are modified using our 
neutron multiplicity correction, good agreement is also 
found between the two results. It is relevant to note that 
there is a significant fraction of the neutrons with energy 
above 10 MeV. 

TABLE V: The muon-induced neutron flux for six sites (in 
units of 10“® cm“^s“^). The total flux is included along with 
those predicted for neutron energies above 1, 10, and 100 
MeV. 


Site 

total 

> 1.0 MeV 

> lOMeU 

> lOOMeV 

WIPP 

34.1 

10.78 

7.51 

1.557 

Soudan 

16.9 

5.84 

4.73 

1.073 

Kamioka 

12.3 

3.82 

3.24 

0.813 

Boulby 

4.86 

1.34 

1.11 

0.277 

Gran Sasso 

2.72 

0.81 

0.73 

0.201 

Sudbury 

0.054 

0.020 

0.018 

0.005 


2. Neutron Production in Common Shielding Materials 

Fast neutrons can also be created by muons passing 
through the materials commonly used to shield a detector 
target from natural radioactivity local to the surrounding 
cavern rock. Fig. d shows the neutron yield in some 
common shielding materials. We have also included a 
simulation for germanium which will prove useful later 
in this paper when we consider the DSR for experiments 
based on this target material. 



FIG. 15: The muon-induced neutron production rate pre¬ 
dicted for some common detector shielding materials. Note 
that minor variations due to neutron back-scattering have 
been neglected in these calculations. 


The fitted functions have the same form as equa¬ 
tion m but with different values for parameters which 
are provided in Table lVll To convert the neutron produc¬ 
tion rate to the total neutron flux, one multiplies equa¬ 
tion ca by the average muon path length which depends 
upon the detector geometry. 


TABLE VI: Summary of the fitting parameters describing the 
muon-induced neutron production rate in common detector 
shielding materials. 


Material 

Po 

Pi 

Lead 

(7.84±2.21) X 10"“ 

0.86±0.05 

Polyethylene (6.89±1.95) x 10"® 

0.86T0.05 

Copper 

(2.97±0.838) X 10"® 

0.87±0.05 

Germanium (3.35di0.95) x 10 ^ 

0.87T0.05 


Generally speaking, muon-induced neutrons produced 
in a detector target or surrounding shield can be actively 
vetoed in coincidence with the primary, depending on the 
veto efficiency and specific detector geometry. Specific 
examples are provided later in this paper. 
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3. Neutron Energy Spectrum 

It is well known that the energy spectrum of neutrons 
produced by muon spallation is uncertain 1^ 1^ 
144 . [43 and that data is scarce, particularly for studies 
relevant for deep underground sites. Nonetheless, with 
reference to Fig. ^1 and our previous discussion of the 
LVD data sample obtained at Gran Sasso, the FLUKA 
simulations reproduce the data well once the data are ap¬ 
propriately corrected for the quenching of proton-recoil 
energy. Recently, Ref. reported a measurement of 
the muon-induced neutron energy spectrum using 190 
GeV/c muon interactions on a graphite target. The neu¬ 
trons were observed by liquid scintillator detectors and 
the neutron energy distribution was determined via time- 
of-flight. The measured angular and energy distributions 
agree well with the FLUKA simulation performed by 
Ref. [ 33 . 



Neutron Energy (GeV) 

FIG. 17: The differential energy spectrum for muon-induced 
neutrons at the various underground sites. The bin width is 
50 MeV. 



FIG. 16: The differential energy spectrum for muon-induced 
neutrons as measured in the LVD experiment before and af¬ 
ter correcting for proton-recoil quenching effects described in 
the text. Following such corrections, the FLUKA simulation 
appears to reproduce the shape of the spectrum well. 


average neutron energy for each site are summarized in 
Table IWl 

TABLE VII: Summary of the fitting parameters describing 
the shape of the differential energy spectrum of muon-induced 
neutrons for the various underground sites. 


Site 

< E > 

ap 

ai a2 

as 

WIPP 

62 MeV 

6.86 

"Tl 2.971 xlO"^'" 

2.456 

Soudan 

76 MeV 

7.333 

2.105 -5.35x10"^® 

2.893 

Kamioka 

79 MeV 

7.55 

2.118 -1.258x10"^'^ 

2.761 

Boulby 

88 MeV 

7.882 

2.212 -2.342x10"^'^ 

2.613 

Gran Sasso 91 MeV 

7.828 

2.23 -7.505x10"^® 

2.831 

Sudbury 

109 MeV 7.774 

2.134 -2.939x10"^® 

2.859 


4- Neutron Angular Distribution 


We derive the neutron energy spectrum for each ex¬ 
perimental site fFig. 11711 from the FLUKA simulation for 
the neutrons produced in the rock and then emerge into 
the experimental hall. The muons are generated locally 
for each site as described in Section II and used as input 
to the FLUKA simulation. 

For each site we provide a convenient parameterization 
based upon the following fitting function: 

AN p-a-oEn 

— = A^{—— + ( 14 ) 

d-tL/n 

where A^ is a normalization constant, ao, ai, a 2 and aa 
are fitted parameters, E„ is the neutron energy, B^(E^) 
is a function of muon energy and E^ is in GeV, 

B^{E^) = 0.324 - 0.641e-° °i4^''. (15) 

This parameterization is consistent with Ref. and 
is valid for E„ >10 MeV. The fit parameters and the 


The angular distribution of neutrons produced in the 
rock by muons is shown in Fig. 1181 As described in 
Refs. dE® , our simulations reproduce the expected for¬ 
ward peak for those neutrons that are produced largely 
through muon spallation whereas the secondary evapo¬ 
ration of neutrons is predominantly distributed isotropi¬ 
cally along the muon track. 

We parameterize the angular distribution according to: 


dN _ _Ae_ 

dcos{0) {1 — cos{6))^o^^A-\-Ce{Efj)' 


(16) 


where Ag is a constant and Be(E^) and Ce(E^) are 
weakly correlated to muon energy and E^ is in GeV. The 
corresponding functions are: 

Be[E^) = 0 . 482 £; 0-°45 


(17) 
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FIG. 18: Simulation of the muon-induced neutron angular 
distribution for neutrons produced relative to the primary 
muon track. 

and 

Ce{E^) = 0.832i;;Oi52. (ig) 

5. Neutron Multiplicity 

The number of neutrons produced per muon interac¬ 
tion is the least known quantity in the production of neu¬ 
trons induced by muons. As discussed previously, the 
average multiplicity in FLUKA is smaller than that of 
the measurements |8lL l4l| . The neutron multiplicity dis¬ 
tributions obtained from our simulations are shown in 
Fig. [T^ 



FIG. 19: Calculated neutron multiplicity at different exper¬ 
imental sites. Shown is a solely FLUKA calculation without 
correcting multiplicity using equation GOJ. 

The parameterization function proposed by Ref. is 
employed: 

^ ( 19 ) 


where Am is normalization constant, M is the multiplic¬ 
ity and is in GeV. We found 

Bm{E^) = 0.321A-°-247, (20) 


C'm(A^) = 318.1e-°°^42ii5^^ (21) 

and 

Dm{E^) = 2 . 026 -°°°®^*^®'®''. ( 22 ) 

The average multiplicity exhibits the expected depen¬ 
dence on muon energy, and thus depth, and is apparent 
in the fit parameters < M > = 3.48, 4.26, 5.17, 6.03, 
6.44 and 7.86 for WIPP, Soudan, Kamioka, Boulby, Gran 
Sasso and Sudbury, respectively. The neutron multiplic¬ 
ity is also dependent on the different target materials. 
We have simulated the neutron distributions using the 
average density and chemical composition appropriate 
for each site and applied a correction to the simulated 
multiplicity according to equation m- The corrected 
multiplicity agrees with the result in Ref. M for Kam- 
LAND. 


6. Neutron Lateral Distribution 

Fig. 1201 shows a FLUKA simulation of the lateral dis¬ 
tribution of neutrons as they emerge from the primary 
muon track in a selection of media. Typically speaking, 
the neutron flux is attenuated by about two orders of 
magnitude at distances larger than 3.5 m from the muon 
track, however, as much as 10% remain at distances as 
large as 2 to 2.5 m. 

IV. MUON-INDUCED COSMOGENIC 
ACTIVITY 

As the cascades of muon-induced reactions propagate 
through detector materials, the production rate of cos- 
mogenic nuclide j at depth A in a detector volume can 
be expressed as: 

i?j(A)=^ni^ J aijkiEk) ■ (j)k{Ek, X)dEk, (23) 

i k 

where n^ is the number of atoms for target element i 
per kg of material in the detector, aijk{Ek) is the cross 
section for the production of nuclide j from the target 
element i by particles of type k with energy E^, and 
(j)k{Ek,X) is the total flux of particles of type k with en¬ 
ergy Efe. The production cross sections are discussed in 
detail in Ref. where the equivalent photon approx¬ 
imation is used. The energy dependence of the corre- 
sponding cross-section can also be described by cr^ (E) = 

oo- EO-^EI. 
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FIG. 20: The fraction of muon-induced neutrons emerging 
as a function of distance from the primary muon track in 
several media. The curves exhibit distinct features relevant 
to neutron production and propagation in the media. At short 
range from the muon track, neutron production increases with 
distance as the nuclear and electromagnetic shower develops, 
however, neutron production weakens after about 50 cm from 
the muon track and propagation/attenuation of neutrons in 
the medium dominates. 

Neutrons can interact with nuclei to produce long-lived 
radioactive isotopes and secondary neutrons. The pro¬ 
duction rate can be estimated as 

R{h) = J $„(T;, h)-Na- an{E)dE, (24) 

where ^n{E,h) is the flux of neutrons on the detector 
at depth h, is the number of atoms of the target 
and cr„ (E) represents the production cross section of the 
neutron reaction The neutron flux <i)„(E) de¬ 

pends on the neutron energy. Note that the long-lived ra¬ 
dioactive isotopes produced by neutrons near the Earth’s 
surface are the dominant product of the muon-induced 
cosmogenic radioactivity. 

The production of cosmogenic radioactivity depends 
strongly on the target and it must be evaluated specifi¬ 
cally for an individual experiment. Nonetheless, the pro¬ 
duction rate is proportional to the muon, or neutron, flux 
and the relevant interaction cross-section. The energy de¬ 
pendence of the total cross-section for all muon-induced 
radio-isotopes in liquid scintillator was evaluated assum¬ 
ing the power law |^ . 

atotiE^.) (X (25) 

where a varies from 0.50 to 0.93 with a weight mean 
value < a > = 0.73±0.10 [^. For a target of N atoms 
and cross-section ao at the Earth’s surface, where the 
average muon energy is about 4 GeV, the muon-induced 
cosmogenic radioactivity (Riso) depends on the differen¬ 
tial muon energy spectrum dN^/dE^ at the experimental 
site: 


As a simplihcation, the production rate is written as a 
function of the average muon energy < E^ > 

R^so = (27) 

where is the total muon flux at the experimental site 
and /3 o. 73 = 0.87 ± 0.03 is the correction factor for the 
averaging of E^ |^. For a given detector target, a sim¬ 
ple scaling relation, or Depth-Sensitivity-Relation (DSR) 
factor F, can be derived, 

p _ Riso (Surface) ^ 4 GeV (f>^(Surface) 

Riso(Underground) ^<E^>' (j) ^(Underground)' 

m 

which describes the reduction in muon-induced activity 
as one moves to deeper and deeper sites. Table Enn 
summarizes this effect. 

TABLE VIII: The scaling factor, F, relevant to the Depth- 
Sensitivity-Relation (DSR) developed for the underground 
sites considered in this work. 


Site 

Depth < Efj,> 

(km.w.e.) (GeV) 

F 

WIPP 

1.585 

184 

2140 

Soudan 

1.95 

212 

4600 

Kamioka 

2.05 

219 

5690 

Boulby 

2.805 

264 

19180 

Gran Sasso 3.1 

278 

29270 

Sudbury 

6.011 

356 

1.67x10® 


V. THE DEPTH-SENSITIVITY-RELATION 

In this section we develop the Depth-Sensitivity- 
Relation (DSR) for the major components of the muon- 
induced background, namely the cosmic-ray muons them¬ 
selves, the induced neutron background, and cosmogenic 
radioactivity. In Fig. [2] we show a global view for the 
DSR where we have arbitrarily normalized the DSR Fac¬ 
tor F at the shallower depth characteristic of the WIPP 
site. Generally speaking, the muon flux and induced ac¬ 
tivity is reduced by about one order of magnitude for 
every increase in depth of 1.5 km.w.e.. 

The curves shown in Fig. [2] are indicative of the rel¬ 
ative muon flux and muon-induced activity that will be 
present for a given laboratory site at its characteristic 
depth. The effect of this activity will depend on the 
specihc details of a given detector geometry, including 
shielding, and the goals of a particular experiment. Gen¬ 
erally speaking, muons that traverse a detector unvetoed 
and muon-induced fast neutrons are the primary concern 
deep underground, while long-lived cosmogenic activity 
is usually dominated by activation of detector materials 
at the surface and prior to construction underground. In 
what follows, we describe simulations for a germanium- 
based detector and apply this to develop the DSR for 
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FIG. 21: Relative attenuation factors for the muon and 
muon-induced activities as a function of overburden. The 
curves are normalized, arbitrarily, to unity for the shallower 
depth defined by the WIPP site at 1.585 km.w.e.. Roughly, an 
attenuation factor of about one order of magnitude is achieved 
for every 1.5 km.w.e. increase in depth. 


specific examples in the search for dark matter and neu¬ 
trinoless double beta decay. Particular attention is paid 
to the sensitivity to muon-induced fast neutrons. 


A. Simulation Set-Up for Germanium Based 
Detectors 

To evaluate the response of neutrons in a detector, a 
Monte Carlo simulation code has been developed to sim¬ 
ulate the neutrons generated in different media where 
we rely on the neutron fluxes and distributions gener¬ 
ated and discussed above. The detector geometry, mate¬ 
rial, and electromagnetic interactions are simulated us¬ 
ing GEANT 3 |H3|- Hadronic interactions are simulated 
using the Nucleon Meson Transport Code, NMTC [h^ . 
while transportation of low energy neutrons is achieved 
using GCALOR F^ist neutrons deposit their en¬ 

ergy via elastic scattering and/or inelastic scattering pro¬ 
cesses. We will demonstrate that the former is the main 
concern for dark matter searches since the elastic scatter¬ 
ing process tends to deposit energy in the low energy re¬ 
gion of interest while the latter dominate the background 
through inelastic scattering process owing to the ensuing 
7 -rays produced above the Q-value for double beta de¬ 
cay. In general, inelastic reactions of fast neutrons leave 
the residual nucleus in a highly excited state which subse¬ 
quently decays via y-cascades to the ground state in typi¬ 
cally three or four steps. The initial intensity distribution 
over a very large number of highly excited levels is col¬ 
lected in the first few excited levels which then decay to 
the ground state. In the simulation, the Hauser-Feshback 
theory is used to calculate inelastic scattering cross 
sections for excitation of a given level depending on the 
properties of the ground state and the excited state. This 


theory was first formed by Hauser and Feshback in the 
1950’s and later modified by Moldauer |^. Since then 
many experiments have verified the theory 

In addition to the geometry associated with the detec¬ 
tor and shielding materials, it is important also to define 
the geometry and dimension of the cavern housing the ex¬ 
periment. For example, we have demonstrated that the 
neutron flux incident on the shielding around a detector 
can vary by factors of about 2-3, depending on the cav¬ 
ern size, due to the back scattering of neutrons from the 
cavern walls. As such we specify a cavern size 30 x 6.5 
x4.5 in our simulations. The effects of lead, polyethy¬ 
lene, copper, and target material on neutron production 
and absorption are also important to the neutron simu¬ 
lation. We have seen a large increase (a factor ^10-20, 
depending on the thickness of lead) in the neutron flux 
due to the additional and efficient production of neutrons 
in lead, an effect that has also been identified in Ref. [^ . 
Consequently, the DSR developed in what follows should 
be understood within the boundary conditions described 
for the specific experiments considered. 

B. DSR for Dark Matter Experiments 

Experiments geared toward the direct detection of dark 
matter such as WIMPs (Weakly Interacting Massive Par¬ 
ticles) rely on detector technologies capable of visible en¬ 
ergy thresholds well below 100 keV in order to observe 
the recoil energy induced via WIMP scattering off the nu¬ 
cleus. In order to have sufficient sensitivity to the feeble 
WIMP cross-section, such detectors must also be con¬ 
structed of materials with extremely low levels of natu¬ 
ral radioactivity and be able to discriminate background 
from ionizing y-rays and electrons that can otherwise fog 
a potential WIMP signal. With this discrimination power 
in hand, it remains to assure that nuclear recoil events 
associated with fast neutrons are kept sufficiently rare 
as they present an ineluctable background in the search 
for WIMPs. To date, the most stringent limits on the 
WIMP-nucleon cross-section 1.6 x cm^) have 
been provided by the CDMS-H experiment operating in 
the Soudan mine and it is the goal of next generation 
experiments to improve this sensitivity by several orders 
of magnitude. 

In order to demonstrate the sensitivity of dark mat¬ 
ter experiments to muon-induced fast neutrons we derive 
the DSR for the CDMS-H detector [^, which consists 
of a tower of four Ge (250 g) and two Si (100 g) detec¬ 
tors surrounded by an average of 0.5 cm of copper, 22.5 
cm of lead and 48.6 cm of polyethylene. A 5-cm-thick 
muon veto detector with efficiency > 99.9% encloses the 
shielding. 

The production rate (R) of nuclear recoil events pro¬ 
duced by fast neutrons can be expressed as: 

R^ih) = (29) 
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where rii is the number of atoms for target element i per 
kg material in the detector, neutron energy 

spectrum (equation 11411 1 at depth <Ji{E) is the neutron 
interaction cross section |4^ with element of 

natural Ge, and Fi{E) is an energy-dependent quenching 
function specific to the element of Ge. 

We generate the muon-induced neutrons at the 
rock/carven boundary using the formalism outlined in 
Section III and propagate them through the GDMS II ge¬ 
ometry described above. Since the muon veto in CDMS 
II has an efficiency greater than 99.9% , we are con¬ 

cerned only for the neutrons produced in the rock. We 
have performed our simulations for CDMS II using two 
different shielding configurations. Shielding configura¬ 
tion 1 is that used in the actual experiment with 0.5 cm 
copper, followed by 8.6 cm polyethylene, 22.5 cm lead, 
and 40 cm of polyethylene as the outer neutron shield. 
In shielding configuration 2, we interchanged the thick 
polyethylene and lead shield positions. In this case we 
found a reduction in background by about a factor of 
two over the CDMS-II shield. This reduction occurs ow¬ 
ing to the additional neutrons produced when neutrons 
from the rock interact in the lead shield. Similar obser¬ 
vations have been made in Ref. and Ref. The 
visible recoil energy spectrum induced by the fast neu¬ 
trons is shown in Fig. EH for three different depths and 
along with the spectrum expected for dark matter as¬ 
suming a cross section Up = 10“^®cm^ and a 100 GeV 
WIMP mass. 



FIG. 22: The predicted event rates for spin-independent 
WIMP-nucleon scattering (dotted-line) in Ge assuming a 
WIMP-nucleon cross-section of ap = 10“^®cm^ and a 100 
GeV WIMP mass. Muon-induced neutron backgrounds are 
also displayed for comparison, indicating the need for greater 
and greater depth as experiments evolve in scale and sensi¬ 
tivity. 

Using these results we determine an event rate of 0.9 
events/kg-year in an energy window of 10 to 100 keV for 
the GDMS-II experiment operating at the Soudan mine. 
This rate is reduced to 0.5 events/kg-year after identify¬ 
ing those neutrons that interact with two or more crys¬ 


tals in the GDMS-II tower. Our prediction is consistent 
with the upper bound of 0.94 ± 0.38 events/kg-year that 
can be deduced from the CDMS II collaboration’s limit 
of 0.05 ± 0.02 neutrons detected during an exposure of 
19.4 kg-day or 34 kg-day [g^. In Ref. 0, Kamat 
also simulated the un-vetoed neutron rate in the CDMS- 
II detector and obtained 0.05 ± 0.02 neutrons for the 
19.4 kg-day exposure, in excellent agreement with our 
prediction. 

We can now derive the DSR appropriate to the CDMS- 
II experiment. As shown in Fig. 1231 the experiment’s 
sensitivity would be limited to ap ~ 10“^‘*cto^ due to 
the muon-induced fast neutron flux at Soudan and that 
depths in excess of 5 km.w.e. will be required to push 
beyond ap ^ 10“"^®cm^, unless the neutron flux can be 
suppressed effectively either by further shielding and/or 
active veto. 
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FIG. 23: The Depth-Sensitivity-Relation (DSR) derived for 
the CDMS-II detector geometry for the two shielding conhgu- 
rations described in the text. The muon-induced background 
is dominated by elastic scattering of neutrons depositing vis¬ 
ible energy in a 10 to 100 keV window. Specific points are 
shown, for example, at the depth of the Soudan mine where 
the CDMS-II detector has been operating. Uncertainties re¬ 
flect those present due to uncertainties in the rock composi¬ 
tion and in generating the muon-induced fast neutron flux. 


It was pointed out in Ref. that the nuclear recoil 
event rate in coincidence with a second energy deposi¬ 
tion not associated with nuclear recoils (electrons, pho¬ 
tons, muons etc.) is a factor of 10 more than the rate 
of isolated nuclear recoil. We cannot confirm this state¬ 
ment (a factor of 10) with the GDMS II geometry for 
the neutrons which are produced in the rock and associ¬ 
ated muons that miss the veto. This is likely due to the 
fact that the GDMS II detectors are se gme nted and much 
smaller than that considered in Ref. |65l . In our simu¬ 
lation, heavy charged particles such as pions, kaons and 
protons with kinetic energy greater than 10 MeV that are 
produced together with neutrons by muons in the rock 
are a factor of 10 less than that of neutrons and they 
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are about 40 times smaller in number than the neutrons 
after passing through the shielding. The electrons and 
bremsstrahlung photons that are produced in the rock 
cannot survive the rock and air. The bremsstrahlung 
photons that are produced in the lead through nuclear 
showers induced by neutrons could generate the recoils 
in the detector. However, this contribution is limited by 
the size of the detectors and most of these events are 
multiple crystal events. In the CDMS-II simulation con¬ 
sidered here as an example, we find that only 10% of the 
recoil energy deposited by fast neutrons are coincident 
with secondary particles that can potentially sum with 
the neutron energy deposited in the detector. 

C. DSR for Double Beta Decay Experiments 

To demonstrate the effect of muon-induced activity in 
the search for neutrinoless double beta decay, we consider 
the geometry proposed for the Majorana project 
where a detector module is made up of 57, 1.05 kg, closely 
packed crystals of germanium enriched to 86 percent in 
^®Ge. While details of the shielding for Majorana are 
under development, we consider an inner-most layer with 
10 cm copper, followed by 40 cm lead and an outer-most 
layer of 10 cm polyethylene. An active muon veto outside 
the passive shield is also assumed but one that is limited 
to 90% efficiency to veto nucleons produced inside the 
shielding. 

Both elastic and inelastic reactions of muon-induced 
neutrons are considered, however, unlike the case for 
dark matter, the dominant source of muon-induced back¬ 
ground for the Majorana geometry results from the high- 
energy cascades that evolve from inelastic neutron scat¬ 
tering on the detector and shielding materials that pro¬ 
duce background in the region of interest of the Q-value 
at 2039 keV. The results of our simulation are shown 
in Fig. 1241 with a breakdown of the main contributions 
summarized in Table na 

TABLE IX: Breakdown of the muon-induced background pre¬ 
dicted for the energy range 2000-2100 keV in a Majorana-like 
experiment operating with an overburden characteristic of the 
Gran Sasso Laboratory. 


Reaction 

Events in the ROI 
(events/keV-kg-year) 

''°Ge{n, n'y) 

0.01 

’^‘^Ge{n, n'y) 

0.002 

Cu(n, n'y) 

0.0019 

n'y) 

0.0035 

Elastic Scattering on Ge 0.0036 

Muon hits 

0.0025 

Others 

0.0024 


Here we have performed our simulations assuming that 
the detector was operated at Gran Sasso depth in order 
to directly compare to previous germanium-based exper¬ 
iments situated there. In this case, we find that the to¬ 




FIG. 24: A simulation of the muon-induced background for a 
Majorana-like experiment operating at an equivalent overbur¬ 
den provided by the Gran Sasso Laboratory. In (a) we show 
the full spectrum with an expanded profile in (b) spanning 
the Region-of-Interest (ROI) around the Q-value for neutri¬ 
noless double beta decay at 2039 keV. The peak at 2023 keV is 
characteristic of that produced via the ’^^Ge{n,n''y) reaction. 

tal muon-induced background is about 0.026 events/keV- 
kg-year for Majorana at the depth of Gran Sasso. The 
dominant contribution (82%) to this background results 
from inelastic neutron scattering processes (Ge(n,n'^), 
Pb{n,n'"f) and Cu{n,n'"f)) on the detector target and 
shielding materials. Others (18%) include stopping-muon 
capture on Ge, neutrons that capture on Ge and on Cu, 
and cosmogenic production in-situ. 

It is interesting to set our simulations within the frame¬ 
work of the Heidelberg-Moscow experiment. Comparing 
our simulation to their background model |67l| . we find 
agreement in the prediction of about 0.003 events/keV- 
kg-year due to events escaping the muon veto, however, 
we believe that the muon-induced neutron background 
and in situ cosmogenic production were missed in their 
simulation. We note that the simulated background in 
ref. [H3 is about 20% lower than that actually mea¬ 
sured in the Heidelberg-Moscow experiment. Interest¬ 
ingly enough, this missing 20%, corresponding to 0.022 
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events/keV-kg-year, is precisely what we have found in 
our simulation. 



FIG. 25: The Depth-Sensitivity-Relation (DSR) derived for 
a Majorana-like experiment showing, specihcally, the results 
from this work assuming the detector is operated at a depth 
equivalent to the Gran Sasso Laboratory. The raw event rate 
in the energy region of interest of 0.026 events/keV-kg-year 
can be reduced by a factor of 7.4 by exploiting the detector 
granularity, pulse-shape discrimination (PSD), and detector 
segmentation. The upper curve displays the background sim¬ 
ulated in the case that no active neutron veto is present and 
the lower curve indicates the reduction that would ensue if an 
active neutron veto were present that is 99% efficient. 

The results of our simulations can be used to derive 
the DSR for Majorana as shown in Fig. |2H1 The neu¬ 
tron induced background can be reduced by about a fac¬ 
tor of 7.4 in Majorana owing to the use of crystal-to- 
crystal coincidences and the use of pulse-shape discrim¬ 
ination and segmentation. Nonetheless, to achieve the 
target sensitivity of next generation double-beta decay 
experiments, 0.00025 events/keV-kg-year corresponding 
to the background level required to reach sensitivity to 
the atmospheric mass scale of 45 meV Majorana neutrino 
mass, the muon-induced background must be reduced by 
roughly another factor of 100. This can be achieved only 
by operating such a detector at depths in excess of 5 
km.w.e., otherwise an active neutron veto would need to 
be implemented with an efficiency in excess of 99%. 

D. (Q:,n) Background 

Once the depth requirement is satisfied, a proper shield 
against (a,n) neutrons from the environment becomes 
necessary. We use the standard rock and the measured 
neutron flux (3.78x 10 “®cto“^s“^ ) at Gran Sasso 

assuming that all underground labs have the same order 
of neutron flux to establish the shielding requirement for 
(a,n) neutrons. This flux corresponds to an average of 
about 2.63 ppm Q.Ti ppm ^^^Th activity in 

Gran Sasso rock and 1.05 ppm q gy ppj^ 232rp]^ 

activity in Gran Sasso concrete |70| . The neutron energy 


spectrum depending on the rock composition is shown in 
Fig. 1^ As can be seen, the total neutron flux is about 
three orders of magnitude higher than that of neutrons 
from the rock due to muon-induced processes, but the 
energy spectrum is much softer. 



FIG. 26: The neutron energy spectrum arising from (Q,n) 
reactions due to radioactivity in the rock. We predict a harder 
energy spectrum in Gran Sasso rock relative to standard rock 
owing to the presence of carbon and magnesium. 

To demonstrate the neutron flux and energy spectrum 
at different boundaries we show the rock/cavern neutron 
flux and energy spectrum with a shielding for Majorana 
described earlier in Fig. for the depth of Gran Sasso 
as an example. 



FIG. 27: The energy spectrum for fast neutrons produced 
by (a, n) reactions in the rock compared to those induced by 
muon interactions in the rock with and without shielding. The 
lower energy neutrons (< 10 MeV) are quickly absorbed using 
polyethylene shielding, however, the high energy portion of 
the muon-induced neutron flux persists. The addition of lead 
shielding adjacent to a detector can also create an additional 
source of muon-induced neutrons. 

Note that the (a, n) neutrons from the rock are quickly 
attenuated to the level of the muon-induced neutrons be- 
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low 10 MeV with rather moderate shielding whereas the 
higher energy muon-induced neutrons are essentially un¬ 
affected. We note also the increase in the muon-induced 
neutrons with the addition of lead shielding owing to the 
additional neutron production in the heavy target. Con¬ 
sequently, the high energy muon-induced neutron back¬ 
ground is the dominant concern given adequate shielding 
for the lower energy {a,n) neutrons. 

We show the shielding requirement for (a, n) neutrons 
as a function of the thickness of polyethylene in Fig. 1281 
in terms of the sensitivity of dark matter and double 
beta decay. Polyethylene shielding 30 to 40 cm thick is 
required for next generation experiments using Ge in the 
search for neutrinoless double beta decay while about 60 
cm is required for dark matter searches. 
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VI. SUMMARY 

We have provided a comprehensive study of the 
cosmic-ray muon flux and salient distributions as a func¬ 
tion of depth and specific to a set of existing underground 
laboratories around the globe. We have applied these dis¬ 
tributions to simulate the induced background at various 
underground sites and, where possible, made direct com¬ 
parison to the available experimental data in order to as¬ 
sess the accuracy of our predictions. A Depth-Sensitivity- 
Relation has been developed and applied to examples of 
germanium-based detectors used in the search for cosmo¬ 
logical dark matter and neutrinoless double beta decay. 

The cosmic-ray muon flux is well described by a simple 
exponential law over a broad range in depth extending 
from about 1 to 8 km.w.e. We have defined depth in 
terms of the total muon flux obtained at an equivalent 
vertical depth to a site with flat overburden. This re¬ 
moves some of the confusion regarding the average depth 
often quoted for laboratories sited beneath mountains 
where the measured total muon flux is ^ 15 to 20% 
greater than what would be predicted based upon the 
average depth alone. 

Good agreement can be found between the output of 
FLUKA simulations and the available experimental data 
on muon-induced fast neutrons provided one accepts our 
argument to correct the LVD data on both flux and en¬ 
ergy distribution due to quenching effects. In that case 
we find that our simulations reproduce the data well, al¬ 
beit with an overall normalization for the total neutron 
flux that appears to be underestimated by ~ 35%. This 
normalization appears to be greatly improved when one 
corrects the output of the FLUKA simulation to agree 
with experimental data on neutron multiplicity. Glearly, 
more data on the fast neutron energy spectrum and mul¬ 
tiplicity induced by muons would be valuable to further 
bench-mark and tune the FLUKA simulations. 

Our example DSR for dark matter searches is de¬ 
veloped based on a model for the GDMS detector and 
demonstrates that depths in excess of ~ 5 km.w.e. will 
be required in order to circumvent background from the 



Thickness (cm) 

FIG. 28: The (Q,n)-induced background versus polyethylene 
shielding thickness for the CDMS-II and Majorana detec¬ 
tors considered in this work. The upper limit on the spin- 
independent WIMP-nucleon cross section obtained by the 
CDMS II Collaboration is shown in the upper panel for 
comparison along with that predicted for the muon-induced 
neutron background at Soudan and Sudbury. The lower 
panel includes our predicted value for the background in the 
Heidelberg-Moscow experiment (KKDC) na before and after 
the reduction obtained using detector granularity, pulse-shape 
discrimination (PSD), and detector segmentation. 


elastic scattering of fast neutrons contaminating the low- 
energy region of interest for recoiling WIMPs. A similar 
conclusion can be made for neutrinoless double beta de¬ 
cay, modeled after the Majorana detector, where back¬ 
ground following nuclear excitation due to the inelastic 
scattering of fast neutrons is the main culprit. Shallower 
depths make such experiments feasible provided the fast 
neutron flux can be adequately shielded and/or actively 
vetoed. The muon and muon-induced activity increases 
by approximately one order of magnitude for every de¬ 
crease in depth of 1.5 km.w.e.. 

The program developed here has been applied to these 
specific types of experiments and detector geometries, 
however, the distributions presented in parameterized 
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form can now be used as input to new simulations and 
background studies in other detectors of interest. The 
program could also be easily extended to underground 
sites under development that have not been considered 
in this work. More recently, we have begun an experi¬ 
mental program to verify some of our specific predictions 
by irradiating a Ge-detector with fast neutrons. Prelim¬ 
inary results indicate that the data agree well with our 
specific predictions for the Majorana detector. The de¬ 
tails of that work is beyond the scope of this paper and 
will be communicated separately. 
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